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ABSTRACT: 


For  modelers  evaluating  the  impact  of  dam,  levee,  and  structure  failures  the  need  is  to  be  able  to 
determine  the  flood  height  and  timing.  A  two-dimensional  (2-D)  model  needs  to  be  able  to  reproduce  this 
flood  wave  along  the  channel  and  over  dry  ground.  This  report  details  the  testing  of  the  2-D  shallow- 
water  module  of  the  ADaptive  Hydraulics  (ADH)  model  for  surges  over  initially  dry  ground.  ADH 
utilizes  an  unstructured  computational  mesh  that  is  automatically  refined.  Other  modules  in  ADH  include 
three-dimensional  (3-D)  Navier  Stokes  (with  and  without  the  hydrostatic  pressure  assumption)  and 
groundwater  flow.  Testing  is  conducted  in  comparison  to  physical  flume  results  for  two  test  cases.  The 
first  test  case  is  for  a  straight  flume  and  the  second  contains  a  reservoir  and  a  horseshoe  channel  section.  It 
is  important  that  the  model  match  the  timing  of  the  surge  as  well  as  the  height.  In  both  cases  the  ADH 
compared  closely  with  the  flume  results. 
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1  Introduction 


Traditionally,  one-dimensional  (1-D)  hydrodynamic  models  have  been  used 
to  predict  the  impact  of  dam  break,  lock  failure,  levee  breach,  etc.  A  1-D  model 
assumes  uniformity  in  velocity  and  water-surface  elevation  over  the  cross 
section.  Since  the  actual  variation  is  how  energy  is  removed  from  the  system,  a 
1-D  model  will  require  estimation  of  energy  loss  coefficients  for  expansion  and 
contraction.  Generally  these  coefficients  cannot  be  found  by  tuning  since  no 
similar  event  in  the  drainage  basin  has  occurred.  The  water-surface  predictions 
are  sensitive  to  these  coefficients.  Two-dimensional  (2-D)  models  represent  the 
lateral  as  well  as  the  longitudinal  variation  in  velocity  and  water  surface  and  so 
are  able  to  reproduce  the  energy  loss  with  much  less  uncertainty. 

The  U.S.  Army  Corps  of  Engineers  (USACE)  needs  a  2-D  model  that  will 
supply  hydrodynamic  variation  laterally  and  longitudinally  in  sections  with 
geometric  features  such  as  bridges,  hydraulic  structures,  and  potential  failures. 
The  model  should  also  be  able  to  reproduce  levee  breach  and  flooding  events 
with  a  minimum  of  coefficient  tuning.  The  purpose  of  this  investigation  is  to 
develop  a  2-D  modeling  system  that  can  be  rapidly  applied  to  a  site  and  supply 
accurate  hydrodynamic  results. 

The  ADaptive  Hydraulics  (ADH)  model  has  features  that  place  it  as  the 
appropriate  code  to  modify  and  use  for  modeling  of  surges  over  initially  dry  land. 
The  primary  feature  that  ADH  contains  that  sets  it  apart  from  other  codes  is  in  its 
ability  to  automatically  refine  the  mesh  as  demanded  by  the  hydrodynamics. 
Since  dam  and  levee  breaches  will  not  generally  allow  the  user  to  know  ahead  of 
time  where  resolution  will  be  needed,  it  is  important  the  model  itself  can  find  and 
refine  these  areas  of  hydrodynamic  sensitivity.  A  second  important  attribute  of 
ADH  is  that  it  can  run  on  many  platforms  from  desktop  computers  to  parallel 
processing  supercomputers.  The  code  was  designed  from  the  beginning  to  be  a 
parallel  code.  If  the  size  of  the  problem  outgrows  the  desktop  computer  it  can  be 
moved  to  more  powerful  computers. 

This  report  begins  with  a  description  of  the  ADH.  This  includes  the  manner 
in  which  wetting  and  drying  is  addressed  and  is  followed  by  the  actual  equations 
that  are  solved.  The  next  section  demonstrates  ADH  results  compared  with 
laboratory  results.  The  final  section  presents  conclusions  from  the  investigation. 
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2  Model  Description 


ADH  is  a  finite  element  code  that  presently  can  simulate  3-D  unsaturated 
groundwater  flow,  3-D  Navier  Stokes  flow  (both  hydrostatic  and  nonhydrostatic), 
and  2-D  shallow- water  flow.  The  code  uses  linear  interpolation  for  all  the 
variables.  ADH  uses  tetrahedral  elements  in  3-D  and  triangle  elements  in  2-D. 
ADH  is  capable  of  mesh  adaptation  so  that  it  can  automatically  split  elements  to 
produce  a  more  accurate  representation  of  the  hydrodynamics.  Later  it  may 
recombine  elements  if  the  hydrodynamics  no  longer  demands  this  resolution. 

This  is  a  tremendous  advantage  for  the  user.  The  modeler  need  only  create  a 
mesh  that  accurately  captures  the  topography  of  the  domain.  The  model  can  then 
refine  the  mesh,  if  necessary,  to  create  an  accurate  flow  field.  Without  this 
feature  the  user  would  have  to  either  create  a  highly  resolved  mesh  everywhere, 
or  the  user  will  have  to  modify  the  mesh  by  trial  as  a  run  progresses. 

ADH  uses  an  inexact  Newton  technique  to  develop  the  linear  set  of  equations 
from  the  finite  element  statement.  These  algebraic  equations  are  solved  for  the 
dependent  variables  of  the  shallow- water  equations  ( h,u,v )  using  the  BiCGSTAB 
algorithm  (see  Kelley  1995)  and  a  set  of  preconditioners.  The  Newton  solver  is 
also  utilized  for  wetting-drying  problems. 

A  mesh  of  smaller  segments,  called  elements,  represents  the  problem  domain. 
The  independent  variables  (x,y)  and  the  dependent  variables  (h,u,v)  are 
interpolated  in  each  of  these  elements.  The  equations  of  motion  and  continuity 
are  integrated  on  combinations  of  these  elements.  Portions  of  the  overall  domain 
are  going  to  be  dry,  and  other  portions  will  be  wet.  One  could  have  elements  be 
either  wet  or  dry,  that  is,  turned  on  or  turned  off.  However,  this  on-off  switch  can 
result  in  model  instability.  It  also  would  require  a  fine  mesh.  Instead,  ADH 
integrates  the  shallow-water  equations  over  the  wet  portion  of  an  element.  A 
single  element  that  contains  the  waterline  will  only  partially  contribute  to  the 
shallow-water  equations.  In  this  way  the  waterline  is  found  indirectly  through  the 
iterative  process.  This  adds  more  nonlinear  behavior  to  the  system  and  is  likely  to 
increase  the  number  of  nonlinear  iterations  required  for  each  time-step.  However, 
it  allows  the  model  to  find  the  waterline  at  the  same  time  that  the  flow  variables 
are  calculated. 

The  following  section  provides  the  shallow-water  equations  as  used  in  ADH. 
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g  =  acceleration  due  to  gravity 
p  =  fluid  density 
zb  =  channel  bed  elevation 
P  =  pressure  at  the  water  surface 
n  =  Manning’s  roughness  coefficient 
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Co  =  dimensional  constant  (C0  =  1  for  SI  units  and  C0  =  1.486  for  English 
units) 

The  Reynolds  stresses  are  determined  using  the  Boussinesq  approach  to  the 
gradient  in  the  mean  currents: 
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where  vt  =  kinematic  eddy  viscosity  (which  varies  spatially). 

The  ADH  shallow- water  equations  are  placed  in  conservative  form  so  that 
mass  balance  and  the  balance  of  momentum  and  pressure  are  identical  across  an 
interface.  This  is  important  in  order  to  match  the  speed  and  height  of  a  surge  or 
hydraulic  jump. 

The  equations  are  represented  in  a  finite  element  approach.  The  quality  of  the 
numerical  solution  depends  on  the  choice  of  basis  or  trial  function  and  the  test 
function.  The  trial  function  determines  how  the  variables  are  represented  and  the 
test  function  determines  the  manner  in  which  the  differential  equation  is 
enforced.  In  the  Galerkin  approach  the  test  functions  are  chosen  to  be  identical 
with  the  trial  functions.  When  the  flow  is  advection-dominated,  the  Galerkin 
approach  produces  oscillatory  behavior.  The  basic  problem  is  that  the  Galerkin 
form  of  the  test  function  cannot  detect  the  presence  of  a  node-to-node  oscillation 
and  so  allows  this  spurious  solution.  The  approach  used  in  ADH  is  to  enrich  the 
standard  Galerkin  test  function  with  an  additional  term  that  can  detect  and  control 
this  spurious  solution. 

The  Petrov-Galerkin  method  used  here  is  based  on  elemental  constants  for 
coefficients.  This  reduces  the  stabilization  to  the  nonconservative  form.  This  is 
not  a  problem  since  the  stabilization  is  only  applied  within  the  elements  and  uses 
the  Galerkin  test  function  to  enforce  “flux”  balance  across  element  edges.  For 
illustration,  consider  the  shallow-water  equations  in  nonconservative  form 
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where,  A  =  and  B  =  .  The  trial  functions  (or  interpolation/basis 

functions)  are  the  Lagrange  polynomials.  These  are  piecewise  linear  functions 
that  are  continuous  across  element  boundaries.  Spatial  derivatives,  however,  are 
not  continuous  across  these  element  edges.  Each  of  the  dependent  and 
independent  variables  is  interpolated  via  these  trial  functions.  For  example, 

N 

u{x)=  Yj  <f>j  ( x)uj ,  means  that  the  approximate  solution,  u(x) ,  is  made  up  of 

M 

the  product  of  the  trial  function  for  node  j  and  the  nodal  value  at  that  location. 
The  test  function  is  chosen  as: 
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/  =  (Qt,)^2 ,  i.e,  the  square  root  of  the  element  area 
V  =  («,v) ,  the  element  average  velocity  components 
A t  =  the  time-step  size 

The  finite  element  statement  becomes: 
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where,  the  subscript  l  indicates  the  finite  element  approximation.  The  Petrov- 
Galerkin  contributions  are  integrated  on  the  interior  of  the  elements,  but  not 
across  element  edges.  This  contribution  stabilizes  the  Galerkin  approach.  This 
scheme  utilizes  a  single  scaling  factor  a.  This  is  different  from  the  scheme 
reported  in  Berger  and  Stockstill  (1995).  That  scheme  involved  scaling  each 
eigenvalue.  However,  this  method  does  not  converge  using  the  iterative  solver  in 
ADH.  Instead,  a  single  value  scaling  (Equation  1 1)  is  used. 
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3  Model  Tests  and  Results 


A  numerical  model  needs  to  be  produced  in  a  careful  manner.  The  underlying 
approach  must  be  conservative  so  that  mass  and  the  momentum  balance  are 
correctly  represented.  If  not,  the  speed  and  location  of  a  shock  will  be  wrong. 
These  comparisons  with  physical  flume  results  will  allow  a  check  on  the 
numerical  code.  The  first  case  is  a  straight  flume  in  which  a  dam  is  rapidly 
removed.  This  produces  a  surge  downstream  over  initially  dry  ground  and  a 
rarefaction  wave  moving  upstream.  The  second  case  also  involves  an 
instantaneous  dam  removal  and  a  surge  moving  over  dry  ground.  This  case 
includes  a  curved  section  of  channel.  The  straight  channel  is  a  1-D  problem  while 
the  curved  channel  represents  a  more  general  2-D  test  case. 


Case  1:  Straight  Flume 

The  first  case  is  a  comparison  with  results  from  a  study  done  at  ERDC’s 
Waterways  Experiment  Station  (U.S.  Army  Engineer  Waterways  Experiment 
Station  1960;  1961).  This  study  is  of  a  dam  break  in  a  straight  flume.  The  flume 
was  121 .9  m  (400  ft)  in  length  and  had  a  width  of  1 .2  m  (4  ft).  The  flume  had  a 
slope  of  0.0015  m  (0.005  ft).  The  dam  is  situated  61  m  (200  ft)  into  the  flume. 
Water  is  pooled  upstream  of  the  dam  prior  to  the  test.  The  initial  water  surface 
upstream  of  the  dam  is  flat  with  the  depth  at  the  dam  of  0.3048  m  (1  ft). 
Downstream  the  bed  is  dry.  The  water  level  is  recorded  at  several  stations  both 
upstream  and  downstream  of  the  dam.  When  the  dam  is  removed,  a  rarefaction 
wave  propagates  upstream  and  a  surge  downstream.  For  this  report,  the  water 
depth  (overland  head)  at  stations  160,  191,  200  (dam),  225,  275,  345  are 
compared  between  ADH  and  the  physical  model  flume.  The  station  numbers 
indicate  the  location  from  the  upstream  edge  of  the  model.  Stations  160  and  191 
are  upstream  of  the  dam;  sta  200  is  at  the  dam  location;  and  stations  225,  275, 
and  345  are  downstream  of  the  dam.  Figure  1  shows  the  initial  water  surface,  the 
bed,  and  channel  stations.1 


1  Figures  1-7  in  this  report  express  units  of  measurement  in  feet.  To  convert  feet  to 
meters,  multiply  the  number  of  feet  by  0.3048. 
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The  numerical  mesh  contains  4,000  elements  and  2,505  nodes.  The 
downstream  model  boundary  extends  30.5  m  (100  ft)  beyond  that  of  the  flume. 
This  allows  the  surge  to  pass  out  of  the  model  without  contaminating  the  results 
at  the  flume  stations.  The  numerical  model  for  this  case  is  treated  as  a  closed 
region  in  which  the  water  stays  within  the  domain  for  the  entire  run.  The 
upstream  and  downstream  boundaries  are  treated  as  impervious  wall.  Otherwise 
the  numerical  flume  matched  that  of  the  physical  flume. 

The  model  parameters  are  given  in  Table  1.  The  Kinematic  Eddy  Viscosity 
(EVS)  card  is  used  to  set  the  kinematic  eddy  viscosity  and  the  Molecular 
Viscosity  (MU)  card  sets  the  molecular  viscosity.  The  Manning’s  Roughness 
Unit  Conversion  (MUC)  card  is  used  for  unit  conversion  with  the  Manning’s 
roughness  formulation.  It  has  a  value  of  1  for  SI  units  and  1.486  for  English 
units.  Density  (RHO)  and  Gravity  (G)  set  the  values  for  density  and  gravity  in  the 
desired  units.  The  Drying  Limits  (DTL)  is  used  for  wetting  and  drying.  The 
Maximum  Levels  of  Refinement  (ML)  and  Shallow-Water  Refinement  Tolerance 
(SRT)  card  are  used  for  adaptation.  The  ML  card  sets  the  number  of  times  the 
grid  can  adapt  and  the  SRT  card  sets  the  error  level  at  which  the  grid  will  adapt. 

A  value  of  zero  indicates  that  the  same  result  was  obtained  with  or  without 
adaptation  on  this  particular  test  case. 


Table  1 

Test  Conditions  for  Straight  Flume  Dam  Break  Case 

Condition 

Values 

Kinematic  Eddy  Viscosity  (EVS) 

0.09  0.09  0.09 

Molecular  Viscosity  (MU) 

0.0001 

Manning’s  Roughness  Unit  Conversion  (MUC) 

1.486 

Density  (RHO) 

1.94 

Gravity  (G) 

32.2 

Drying  Limits  (DTL) 

0.01  0.04 

Maximum  Levels  of  Refinement  (ML) 

0 

Shallow-Water  Refinement  Tolerance  (SRT) 

100 

Time-step 

0.25 

Manning’s  n 

0.009 

The  time-step  is  determined  by  decreasing  the  size  until  the  results  no  longer 
changed.  The  results  are  considered  converged  with  a  time-step  of  0.25  sec.  The 
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physical  model  study  suggests  a  Manning’s  n  value  of  0.009,  which  is  used  in  the 
numerical  model. 

Figures  2  through  7  show  the  test  results  for  stations  160,  191,  200,  225,  275, 
and  345.  These  graphs  show  a  time-history  of  the  water  level  for  both  the 
numerical  model  and  the  flume.  The  arrival  time  of  the  wave  in  ADH  agrees  with 
the  arrival  time  from  the  flume  at  all  stations.  Initially  the  water  level  in  the 
model  seems  to  fall  at  the  same  rate  as  in  the  flume  for  both  the  upstream 
stations.  Near  the  end  of  the  run,  the  water  level  in  the  flume  drops  off  more 
rapidly  than  predicted  in  the  model  (see  Figures  2  and  3).  The  flume  results  at  the 
site  of  the  dam,  sta  200,  show  a  surge  superimposed  on  the  overall  water-surface 
drop.  The  numerical  model  captures  the  overall  drop  but  not  the  superimposed 
surge.  The  maximum  difference  at  this  station  occurs  at  a  time  of  50  sec. 

The  three  downstream  stations  (225,  275,  and  345),  shown  in  Figures  5,  6, 
and  7,  also  have  some  points  that  differ  from  the  physical  model  by  as  much  as 
0.04.  Overall  the  water  level  predicted  by  the  numerical  model  is  in  close 
agreement  with  the  data  from  the  flume.  At  sta  225  (Figure  5)  the  time  of  arrival 
of  the  water-surface  peak  predicated  in  the  numerical  model  agrees  with  that  of 
the  flume  although  there  is  a  3  percent  error  in  the  maximum  water  level.  The 
time  predicted  for  maximum  water  level  at  sta  275  (Figure  6)  is  5  sec  earlier  than 
seen  in  the  flume.  This  is  about  a  14  percent  error.  The  peak  elevation  at  this 
station  also  differs  by  1 1  percent  from  the  flume.  The  last  station,  345  (Figure  7), 
has  only  a  5  percent  error  in  the  maximum  water  level  and  the  predicted  time 
agrees  with  the  flume.  The  rise  in  the  water  level  in  the  numerical  model  after 
125  sec  is  due  to  the  water  collecting  above  the  lower  boundary  of  the  numerical 
model.  This  is  after  the  time  of  comparison  to  the  physical  flume. 


Figure  2.  Depth  comparison  for  ADH  and  flume  experiment,  sta  1 60 
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Figure  4.  Depth  comparison  for  ADH  and  flume  experiment,  sta  200 
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Figure  6.  Depth  comparison  for  ADH  and  flume  experiment,  sta  275 
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Figure  7.  Depth  comparison  for  ADH  and  flume  experiment,  sta  345 


Case  2:  Horseshoe  Curve  Flume 

The  second  case  is  a  comparison  with  a  study  done  by  Bell  et  al.  (1992). 
Figure  8  shows  a  plan  view  of  the  flume  used  in  the  study.  The  flume  is 
composed  of  a  reservoir  section  and  a  channel  section.  The  reservoir  has  a 
bottom  elevation  of  0.089  m  (0.291  ft)  lower  than  the  channel  bed.  The  initial 
water  depth  in  the  reservoir  is  0.305  m  (1  ft)  and  the  channel  is  dry.  The  water  is 
pooled  in  the  reservoir.  The  water  level  in  the  flume  and  the  time  was  recorded  at 
five  of  the  eight  stations.  No  data  were  provided  at  stations  3,  5,  and  7  from  the 
flume  study.  Sta  1  is  located  in  the  straight  channel  section  0.3  m  (1  ft) 
downstream  from  the  dam.  Stations  2,  4,  and  6  are  in  the  curved  section  of  the 
flume.  Sta  2  is  3  deg  into  the  curve;  sta  4  is  90  deg  into  the  curve;  and  sta  6  is 
177  deg  into  the  curve.  Sta  8  is  located  2.29  m  (7.51  ft)  from  the  end  of  the 
curve.  The  numerical  model  parameters  are  given  in  Table  2. 


Table  2  | 

Test  Conditions  for  Horseshoe  Curve  Flume  Dam  Break  Case 

Condition 

Values 

EVS 

0.008  0.008  0.008 

MU 

0.00001 

MUC 

1 

RHO 

1000 

G 

9.81 

DTL 

0.004  0.01 

ML 

2  1 

SRT 

3 

Time-step 

0.03 

||  Manning’s  n 

0.009  | 
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The  grid  used  in  the  numerical  model  is  shown  in  Figure  9.  This  initial  grid 
has  899  elements  and  532  nodes.  ADH  automatically  refines  the  mesh  based 
upon  the  hydrodynamic  error  measured  during  the  simulation  so  the  initial  grid  is 
the  smallest  number  of  elements  and  nodes  for  the  model  test.  The  Bell  et  al. 
(1992)  article  suggests  a  Manning’s  n  value  of  0.0165  based  on  their  numerical 
model  study.  However,  prior  experience  suggests  this  is  far  too  high  for  a 
Plexiglas  flume.  Instead  a  Manning’s  n  value  of  0.009  is  used,  which  is 
consistent  with  other  studies  done  in  smooth  flumes,  e.g.,  Berger  (1993);  Berger 
and  Stockstill  (1995);  and  Stockstill  et  al.  (1997);  and  the  straight  channel  test  in 
this  report.  A  Manning’s  n  value  of  0.009  is  in  line  with  the  recommended  values 
found  in  Chow  (1959)  and  Henderson  (1966). 

The  test  result  comparisons  are  shown  in  Figures  10  through  14.  In  these 
graphs,  symbols  indicate  the  experimental  data  and  a  solid  line  indicates  the 
numerical  data.  The  outer  wave  is  shown  in  blue  and  the  inner  wave  in  red.  At 
stations  1,  6,  and  8  (Figures  10,  13,  and  14,  respectively)  the  predicted  arrival 
time  from  the  numerical  model  agrees  with  the  time  observed  in  the  flume.  The 
predicted  arrival  time  at  sta  2  (Figure  1 1)  is  approximately  one  tenth  of  a  second 
early,  which  is  about  a  1  percent  error.  The  predicted  arrival  time  for  the  inner 
wave  at  sta  4  (Figure  12)  is  early  by  three  tenths  of  a  second  (a  4  percent  error) 
and  the  outer  wave  is  early  by  two  tenths  of  a  second  (a  6  percent  error).  All  of 
these  arrival  times  are  likely  within  the  test  uncertainty  in  the  physical  flume.  The 
predicted  water  level  rises  at  the  same  rate  seen  in  the  flume  at  stations  1,  2,  and 
8.  At  sta  4  the  predicted  water  level  rises  at  a  rate  that  is  a  bit  slower  than  seen  in 
the  flume  for  the  outer  wall,  but  along  the  inner  wall  the  water  level  increases  at 
the  same  rate.  The  difference  between  the  inside  and  the  outside  wall  water- 
surface  elevations  is  small.  At  sta  6  the  water  level  along  the  outer  wall  rises  at 
the  same  rate,  but  along  the  inner  wall  the  water  level  rises  more  slowly  than  seen 
in  the  flume.  The  predicted  maximum  height  agrees  with  the  water  levels  seen  in 
the  flume  at  all  stations  except  sta  2.  In  the  results  from  the  flume,  the  maximum 
water-level  record  at  this  location  is  around  0.043  m  (0.141  ft).  At  the  time  these 
maximum  heights  are  observed  in  the  flume,  the  numerical  model  predicts  a 
height  of  0.067  m  (0.219  ft)  and  by  the  end  of  the  run  the  numerical  model  has  a 
height  of  0.09  m  (0.29  ft). 
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Figure  1 1 .  Depth  comparison  for  ADH  and  flume  experiment,  sta  2 
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Figure  13.  Depth  comparison  for  ADH  and  flume  experiment,  sta  6 
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Figure  14.  Depth  comparison  for  ADH  and  flume  experiment,  sta  8 
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4  Conclusions 


One  of  the  key  characteristics  for  this  type  of  problem  that  ADH  must 
reproduce  is  the  timing  of  the  surge.  If  the  formulation  is  not  fully  conservative, 
the  speed  of  the  surge  will  not  match.  There  are  six  comparison  locations  in  test  1 
and  eight  locations  in  test  2.  Overall  the  two  test  cases  show  ADH  matches  the 
surge  speed  very  well  at  all  stations. 

The  other  key  characteristic  is  the  height  of  the  surge.  Of  the  six  data 
locations  in  test  1,  five  of  the  six  closely  match.  In  the  eight  data  locations  in  test 
2,  six  compare  closely.  ADH  matches  surge  heights  well. 

Parameters  are  selected  based  upon  reasonably  well  known  values.  The 
model  comparisons  are  accomplished  without  need  for  tuning.  This  investigation 
shows  that  ADH  compares  well  with  tests  of  dam  break  and  surge  propagation. 
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